#include "randomc.h"

void TRandomMersenne::RandomInit( uint32 seed )
{
    // re-seed generator
    mt[0] = seed;
    for ( mti = 1; mti < MERS_N; mti++ )
    {
        mt[mti] = ( 1812433253UL * ( mt[mti - 1] ^ ( mt[mti - 1] >> 30 ) ) + mti );
    }

    // detect computer architecture
    union
    {
        double f;
        uint32 i[2];
    } convert;
    convert.f = 1.0;
    // Note: Old versions of the Gnu g++ compiler may make an error here,
    // compile with the option  -fenum-int-equiv  to fix the problem
    if ( convert.i[1] == 0x3FF00000 ) Architecture = LITTLE_ENDIAN1;
    else if ( convert.i[0] == 0x3FF00000 ) Architecture = BIG_ENDIAN1;
    else Architecture = NONIEEE;
}


void TRandomMersenne::RandomInitByArray( uint32 seeds[], int length )
{
    // seed by more than 32 bits
    int i, j, k;
    RandomInit( 19650218UL );
    if ( length <= 0 ) return;
    i = 1;
    j = 0;
    k = ( MERS_N > length ? MERS_N : length );
    for ( ; k; k-- )
    {
        mt[i] = ( mt[i] ^ ( ( mt[i - 1] ^ ( mt[i - 1] >> 30 ) ) * 1664525UL ) ) + seeds[j] + j;
        i++;
        j++;
        if ( i >= MERS_N )
        {
            mt[0] = mt[MERS_N - 1];
            i = 1;
        }
        if ( j >= length ) j = 0;
    }
    for ( k = MERS_N - 1; k; k-- )
    {
        mt[i] = ( mt[i] ^ ( ( mt[i - 1] ^ ( mt[i - 1] >> 30 ) ) * 1566083941UL ) ) - i;
        if ( ++i >= MERS_N )
        {
            mt[0] = mt[MERS_N - 1];
            i = 1;
        }
    }
    mt[0] = 0x80000000UL;
} // MSB is 1; assuring non-zero initial array


uint32 TRandomMersenne::BRandom()
{
    // generate 32 random bits
    uint32 y;

    if ( mti >= MERS_N )
    {
        // generate MERS_N words at one time
        const uint32 LOWER_MASK = ( 1LU << MERS_R ) - 1;       // lower MERS_R bits
        const uint32 UPPER_MASK = 0xFFFFFFFF << MERS_R;        // upper (32 - MERS_R) bits
        static const uint32 mag01[2] = {0, MERS_A};

        int kk;
        for ( kk = 0; kk < MERS_N - MERS_M; kk++ )
        {
            y = ( mt[kk] & UPPER_MASK ) | ( mt[kk + 1] & LOWER_MASK );
            mt[kk] = mt[kk + MERS_M] ^ ( y >> 1 ) ^ mag01[y & 1];
        }

        for ( ; kk < MERS_N - 1; kk++ )
        {
            y = ( mt[kk] & UPPER_MASK ) | ( mt[kk + 1] & LOWER_MASK );
            mt[kk] = mt[kk + ( MERS_M - MERS_N )] ^ ( y >> 1 ) ^ mag01[y & 1];
        }

        y = ( mt[MERS_N - 1] & UPPER_MASK ) | ( mt[0] & LOWER_MASK );
        mt[MERS_N - 1] = mt[MERS_M - 1] ^ ( y >> 1 ) ^ mag01[y & 1];
        mti = 0;
    }

    y = mt[mti++];

    // Tempering (May be omitted):
    y ^=  y >> MERS_U;
    y ^= ( y << MERS_S ) & MERS_B;
    y ^= ( y << MERS_T ) & MERS_C;
    y ^=  y >> MERS_L;
    return y;
}


double TRandomMersenne::Random()
{
    // output random float number in the interval 0 <= x < 1
    union
    {
        double f;
        uint32 i[2];
    } convert;
    uint32 r = BRandom(); // get 32 random bits
    // The fastest way to convert random bits to floating point is as follows:
    // Set the binary exponent of a floating point number to 1+bias and set
    // the mantissa to random bits. This will give a random number in the
    // interval [1,2). Then subtract 1.0 to get a random number in the interval
    // [0,1). This procedure requires that we know how floating point numbers
    // are stored. The storing method is tested in function RandomInit and saved
    // in the variable Architecture. The following switch statement can be
    // omitted if the architecture is known. (A PC running Windows or Linux uses
    // LITTLE_ENDIAN1 architecture):
    switch ( Architecture )
    {
    case LITTLE_ENDIAN1:
        convert.i[0] =  r << 20;
        convert.i[1] = ( r >> 12 ) | 0x3FF00000;
        return convert.f - 1.0;
    case BIG_ENDIAN1:
        convert.i[1] =  r << 20;
        convert.i[0] = ( r >> 12 ) | 0x3FF00000;
        return convert.f - 1.0;
    case NONIEEE:
    default:
        ;
    }
    // This somewhat slower method works for all architectures, including
    // non-IEEE floating point representation:
    return ( double )r * ( 1. / ( ( double )( uint32 )( -1L ) + 1. ) );
}


int TRandomMersenne::IRandom( int min, int max )
{
    // output random integer in the interval min <= x <= max
    int r;
    r = int( ( max - min + 1 ) * Random() ) + min; // multiply interval with random and truncate
    if ( r > max ) r = max;
    if ( max < min ) return 0x80000000;
    return r;
}
